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This Supplementary Information provides details about the spectral extraction (crowding issues 
and the removal of nebular lines), a discussion about the distance to M33, a model for the 0-star 
wind and the measurement of the true photospheric X-ray eclipse width 0, and details about ellip- 
soidal modelling. It also contains five related figures, one related table, and additional references. 
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1 Supplementary methods 



Crowding issues: M33 X-7 lies in a moderately crowded OB association known as HS 13 (ref. 
23). Fortunately, HS 13 has been observed using the Hubble Space Telescope and the WFPC2. 
In ref. 5 it is shown that apart from a close pair (0.2 arcsecond separation) of stars 0.9 arcseconds 
to the southwest, the optical counterpart of M33 X-7 is relatively isolated. For the spectroscopic 
observations from the Gemini-North Telescope, the 0.5 arcsecond wide slit was placed at a position 
angle of 215.6 degrees, which is the angle defined by M33 X-7 and the close pair of neighbour stars, 
which are the only potential source of contamination. 

The detector in the GMOS instrument consists of three 2048 x 4068 EEV CCD chips in a row 
with « 0.5mm gaps in between. The dispersion axis runs along the longer dimension of the mosaic. 
Supplementary Figure \T\ shows spatial profiles of M33 X-7 and the neighbour stars near H/3 for 
three observations that represent relatively "poor" seeing (about 0.8 arcseconds), "average seeing" 
(about 0.6 arcseconds), and "good" seeing (about 0.4 arcseconds). The majority of the spectra were 
similar to the "average" case. Gaussian fits to the spatial profiles and simple numerical integration 
show that the area in the overlap region of the two profiles is typically 5-10% and at most 20% of 
the area under the M33 X-7 profile. For the spectral extraction with the GMOS pipeline, the lower 
extraction aperture was placed slightly to the right of the "dip" between the two profiles, which 
was sufficient to exclude light from the neighbour stars. 



Removing nebular lines: As noted above, M33 X-7 lies in an OB association, and light from the 
surrounding HII region is evident in both the direct images and in the two-dimensional spectra. 
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Unfortunately, the background nebular light varies over small spatial scales and is not removed by 
the regular GMOS spectral extraction routines, which assume that the background light is uniform 
over the entire spatial extent of the two-dimensional spectra. Also, the presence of the close pair 
of neighbour stars makes it harder to use higher order polynomials when fitting the background 
region. 

The SPECRES package in IRAF provides tools to extract point source spectra with complex 
backgrounds 22 . The spatial profile of a point source is distinguished from the quite different profile 
of an extended background source using a Gaussian smoothing kernel along with a Richardson- 
Lucy type of iterative restoration algorithm that performs a maximum likelihood estimation. For 
optimal results, one must precisely specify the positions of the point sources, the variation of the 
spatial position of the point sources with wavelength, and the spatial "Slit Spread Function" (SSF) 
with wavelength. A list of the relative positions of the point sources along the slit was made from 
the Gemini direct images and also the archival HST images. Given the list of relative positions, 
the positions of the sources on a given observation can be found once the position of a reference 
source is measured. The variation of the position of the spectra with wavelength (the "slope") is 
easily measured using the extraction tools in the GMOS package. Unfortunately, the derivation of 
the SSF is somewhat more involved since it depends on the properties of the instrument and on 
the observing conditions. The "specpsf" task in SPECRES was employed to make SSFs from the 
direct images that were taken immediately following each spectroscopic observation and by using 
analytic functions (Gaussians and Lorentzians) with widths that varied with wavelength. 
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Once the point source positions, the slope, and SSFs are determined, one must decide on 
a width for the smoothing kernel that is used to distinguish between point sources and extended 
background sources. There is a trade-off between the width of the kernel and the signal-to-noise 
(SN) in the extracted spectra. A broad smoothing kernel will result in spectra with higher SN, but 
the background subtraction will be poor if the background varies on small spatial scales. On the 
other hand, a narrow smoothing kernel will result in spectra with good background subtraction, but 
at the expense of lower SN in the extracted spectrum. 

Since the observing conditions change from observation to observation, the optimal set of 
parameters for one spectrum may not be appropriate for another spectrum. To optimise the spectral 
extraction for each observation, we wrote scripts that performed the spectral extraction with a wide 
range of variations of the position of the reference source, the spectral slope, the variation of the 
width of the SSF with wavelength, and the size of the smoothing kernel. For each observation, 372 
extractions were performed, and the optimal extraction was determined by cross correlation with a 
synthetic spectrum as the template. 

With the exception of the stronger [O III] line and H/5, the nebular lines were cleanly removed 
in all of the spectra. Supplementary Figure [2] shows a spectrum extracted with the GMOS pipeline 
software and the spectrum extracted with SPECRES. The nebular lines are mostly gone in the latter 
spectrum, although note the higher noise level. 

Since the H/3 emission line was never cleanly removed, we excluded this feature when mod- 
elling the spectrum and computing the radial velocities. As a check on the results, we note that the 
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radial velocity curve derived from the SPECRES extracted spectra (using several Balmer lines and 
He I lines in the cross correlation region 4000-4375 and 4450-4565 A) is virtually the same as the 
radial velocity curve derived from the GMOS pipeline extracted spectra (using the two strong He 
II lines noted in the main text). For example, for the former, we find K 2 = 112.0 ± 7.6 km s _1 and 
for the latter we find K 2 = 108.9 ± 6.4, which are in close agreement with our adopted value in 
Table 2, main text. 

The distance to M33: We estimate the distance to M33 using four indicators. Because of the dif- 
ferential nature of distance measurements, we express our intermediate results as relative distance 
moduli with respect to the Large Magellanic Cloud (LMC), and then place the final result on an 
absolute scale. 

Tip of the Red Giant Branch: Tip of the red giant branch (TRGB) magnitudes for ten fields 
in M33 are given in ref. 26, and a global TRGB value for the LMC is derived in ref. 27. Both 
measurements were carried out using the same calibration of the technique, which includes a small 
(-0.03 mag) metalicity correction. We find A//(M33 - LMC) = 6.22 ± 0.10 mag. 

Cepheids: V and I photometry of 61 Cepheids in M33 with periods from 15 to 80 days 
is presented in ref. 28, which were fitted with LMC period-luminosity (P-L) relations from the 
OGLE project 29 . We apply an updated correction for the metalicity dependence of the Cepheid 
P-L relation, based on the results in refs. 30 and 31, that amounts to +0.06 mag. Thus, we obtain 

A^(M33 - LMC) = 6.18 ± 0.03 mag. 
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RR Lyme: Sixty four RR Lyrae variables in M33 discussed in ref. 32 were used to derive 
an average extinction-corrected V magnitude for the sample of (V) = 25.34 ± 0.07 mag. New 
photometry and spectroscopy for over 100 objects in the bar of the LMC is presented in ref. 33, 
from which an average value of (V) = 19.07 ± 0.06 mag was derived. While these values are 
already corrected for the metalicity dependence of the RR Lyrae, the calibrations used by the 
authors are slightly different (e.g. ref. 34 vs. ref. 35). This results in an additional —0.02 mag 
correction to the results in ref. 32, bringing the relative distance modulus to Ayu(M33 — LMC) = 
6.25 ± 0.09 mag. 

Red Clump: An extinction-corrected /-band magnitude of the red clump in the LMC of 
I = 18.12 ± 0.06 mag is given in ref. 33. An equivalent value for M33 of I = 24.43 ± 0.04 mag 
was determined in ref. 26. Both references use two sets of calibrations for the red clump magnitude 
as a function of metalicity 36 - 37 which results in a small (+0.03 mag) metalicity correction and a 
relative distance modulus of A//(M33 - LMC) = 6.34 ± 0.07 mag. 

The average distance modulus of these four techniques, weighted by their uncertainties, is 
A^(M33 - LMC) = 6.21 ± 0.03 mag. 

We calculate the true distance modulus of the LMC based on the weighted average of three 
recent and independent estimates. First, the study of several detached eclipsing binary systems 
(DEBs) 38 yields an LMC distance modulus of 18.42 ± 0.06 mag. Second, parallaxes of Galactic 
Cepheids obtained with the Hubble Space Telescope yield an absolute calibration of the Cepheid 
P-L relation that, when applied to the LMC, gives a distance modulus of 18.40 ± 0.05 mag (after 
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a correction for metalicity dependence). Lastly, the relative Cepheid distance modulus between 
NGC 4258 and the LMC 30 , coupled with the maser distance to the former 39 , gives 18.41±0.16 mag. 
The average of these values (weighted by their uncertainties) is //(LMC) = 18.41 ± 0.04 mag. 

In conclusion, the combination of A/x(M33- LMC) and jt(LMC) yields jz(M33) = 24.62 ± 
0.05 mag, which corresponds to a distance of d = 840 ± 20 kpc. 

We note that our result is in disagreement with the recent distance determination to M33 in 
ref. 40, who obtain //(M33) = 24.92 ± 0.12 mag, or d = 960 ± 50 kpc based on a DEB. This is 
a puzzling result, because it implies a relative distance modulus between the LMC and M33 based 
on DEBs of 6.50 ± 0.13 mag. The DEB technique has been applied to a similar system in M31 
(ref. 41) and in that case, the relative distance modulus between M31 and the LMC is in excellent 
agreement with other methods. 

Modelling the O-star wind and measuring the photospheric X-ray eclipse width in M33 X-7: 

The mean X-ray luminosity of M33 X-7 is 8 x 10 37 erg s _1 . For an inclination of 75°, and using a 
limb darkening law 42 given by I = Iq(0.5 + 0.75 cosi), the isotropic luminosity is 2.2 x 10 38 erg 
s _1 . For a nominal 10% accretion efficiency, this corresponds to a mass accretion rate of 1.7 x 10 18 
g s _1 = 2.7 x 10 -8 M Q yr _1 . We use this estimate to constrain the O-star wind. 

The X-ray intensity following egress increases steeply by roughly a factor of 15 (~ 0.01 
counts s _1 versus « 0.15 counts s _1 ; Supplementary Figure [3]). Using the out-of-eclipse spectral 
models in ref. 5 and a metalicity of 10% of the solar value, the column density of gas required to 
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attenuate the X-ray flux by a factor of 15 is xs 0.7 g cm 2 . 

As a function of radius r, the velocity v(r) of a radiatively driven O-star wind has a profile 
of the form 43 

v(r) = Voo (l - , (1) 

where is the asymptotic speed of the wind at large radius, R 2 is the radius of the O-star, and (3 is 
an index which typically has a value in the range 0.8 to 1.2. In the following, we take = 2000 
km s -1 , a typical value for O-stars, (3=1, and the value of R 2 given in Table 2 (main text): 
19.6 R e = 1.36 x 10 12 cm. With a = 42 A R Q = 2.95 x 10 12 cm (Table 2, main text), the wind 
velocity at the radius corresponding to the location of the black hole (r = a) is then v(a) = 1081 
km s -1 . 

The mass accretion rate onto the black hole may be roughly estimated via the Bondi-Hoyle 
formula (e.g. ref. 44), 

2 

= 1.75 x 10 32 p(a) g s _1 , (2) 

where p(a) is the density of the wind at r = a. This formula is only approximately correct 
and we should include an unknown coefficient to allow for uncertainties. However, numerical 
simulations 45 show that the coefficient is not very different from unity. 

Equating the expression for M in equation (2) to the previously estimated mass accretion 
rate of 1.7 x 10 18 g s -1 , we find p{a) 1.4 x 10~ 14 g cm~ 3 . The mass loss rate in the O-star wind 
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is then 

M w « 4nap{a)v{a) = 2.6 x 10~ 6 M yr" 1 . (3) 

For comparison, the mass loss rate derived from line-force computations 46 is log M w = —6.274 ± 
0.226 (5.32 x 1CT 7 M yr -1 ), assuming an O-star luminosity and mass given in Table 2 (main 
text), a temperature between 34,000 and 36,000 K, a metalicity between 0. 1 and 0.3 times the solar 
value, and 2.0 < v^/v^c < 2.6. This mass loss rate is a factor of 5 smaller than the value we 
derived. However, the plots that compare the predicted mass loss rates with observed values for 
LMC and SMC stars (Figs. 6 and 7 in ref. 46) have a scatter of up to 0.6 dex (e.g. nearly a factor of 
4). We therefore do not consider the difference between the mass loss rate we derived and the one 
derived from the line-force computations significant. Finally, using M w = 2.6 x lO~ 6 M yr _1 , 
the wind density at radius r is 

/ \ M w 6.6 x 10 10 _ 3 
4nr 2 v(r) ryr — R2) 

We have assumed a spherically symmetric wind, which is perhaps not a very good approximation 47 , 
but it is hard to know how to improve on this. 

The photosphere of the O-star is located at a radius R p such that the radial column density 
from infinity down to R p is equal to 0.54 g cm~ 2 , which is the column depth to the r = 1 surface 
in the best-fitting model atmosphere: 

/ p( r )dr = 0.54 gem' 2 . (5) 

J Rp 

Solving, we find that R p is extremely close to the O-star radius of 19.6 R Q . Thus, the presence of 
the wind has a negligible effect on the radius of the photosphere. 



The X-ray eclipse light curve, however, is strongly affected by the wind. Let us define the 
critical radius of the eclipse R e to be the impact parameter at which the line-of-sight column density 
is equal to 0.7 g cm~ 2 . Thus 

J°° z p (JRI + z 2 ) dz = 0.7 g cm" 2 , (6) 

where Z x = J a 2 — R\ gives the position of the X-ray source. Solving this condition numerically, 
we obtain an eclipse width of = 51°. 

There are uncertainties in many of the parameters we have used. Let us allow for all the 
uncertainties by varying M w by a factor of 1.5 either way around the estimate in equation (3). 
Over this entire range of M w , R p continues to remain very close to R 2 , so the wind has no effect 
on the photosphere. At the lower end of the range, we obtain for the eclipse radius = 49°, while 
at the upper end, we find = 55°. The observed width of = 53 ± 2.2° is bracketed by these 
values. 

According to this model, the photosphere subtends an angle of 46.3° in the light curve, while 
the eclipse, which is defined by the positions at which the X-ray flux falls by a factor of 15, subtends 
an angle in the range 49° — 55°. Thus, the ratio of the eclipse angle to the photospheric angle lies 
in the range 1.06 to 1.19. Turning this around, since the observed eclipse angle is 53°, we estimate 
the photospheric angle to lie in the range 44° — 50°. 

The X-ray light curve itself suggests an angle close to 46° (Supplementary Figure [3]). We 
fitted two line segments to the ACIS count rate data and found that the transition from the flat 
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bottom to egress occurs between = 40° (±1°) and = 46° (±1°), depending on how the data 
are binned and how they are fitted (i.e., x 2 versus median fitting). Thus, in modelling the binary 
system, we set = 46 ± 1° (Table 2, main text). Although we are confident that = 46° is the 
appropriate value to use, we nevertheless show below that derived geometrical parameters change 
only modestly for values of between 40° and 50°. 

Modelling and error estimation: We used the x 2 statistic to evaluate the goodness-of-fit between 
the observed light and velocity curves and their model counterparts: 

2 _ ^fi y(x i ;a 1 ...a w ) - yf | ™ yjxj] ai...a 10 ) - yf 
Xlight — 2^ -.2 + 2^ -.2 

i=l u i i=l u i 
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Here, the notation yf yf yf yf are the observed B, V, g', and r' magnitudes at orbital phase Xi, 
respectively, a { is the uncertainty on the measurement at x { , and y(xi, ai...a w ) is the model value 
at Xi. We also have three observed quantities that constrain the available parameter space, namely 
the duration of the X-ray eclipse 0, the radius of the secondary star R 2 , and the projected rotational 
velocity of the secondary star V rot smi. For each model that is computed, these three quantities 
can be determined. Hence, there are three additional contributions to x 2 ' 
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Our total x 2 is then 

Xtotal = Xlight + XrV + Xe + Xr 2 + Xrof (12) 

In any optimisation procedure, there is always the issue of assigning relative weights to different 
data sets. After our extensive initial searches of parameter space resulted in a good solution, we 
scaled the uncertainties on each data set (e.g., four light curves and one radial velocity curve) so 
that the total x 2 of the fit was equal to N — 1 for each data set separately. The required scalings 
and number of observations for each data set were 0.8715 for 30 observations in B, 0.9505 for 70 
observations in V, 1 .0064 for 39 observations in g', 1 .0449 for 24 observations in r', and 0.9452 for 
22 radial velocity observations. After the scaling, the optimiser codes were run again to produce 
the final set of parameters. The resultant x 2 of the best-fitting solution is ~ 180. 

We computed the confidence limits of the fitted and derived parameters using a brute force 
method. When a model was computed, we saved the value of fitted parameters, the derived pa- 
rameters (e.g., the mass of the compact object, surface gravity of the companion, etc.), and the 
X 2 of the fit. After a suitably large number of runs of the genetic optimiser and the grid search 
optimiser, we computed models and x 2 values almost everywhere near the global x 2 minimum in 
the 10-dimensional parameter space. It was then a simple matter to project one dimension of the 
X 2 hypersurface along any parameter of interest. Supplementary Figure |4] shows x 2 vs. parameter 
value for twelve fitted and derived parameters. A sixth order polynomial was fitted to each curve 
to determine the value of the parameter at the minimum x 2 , and to compute the formal 1, 2, and 3a 
confidence intervals (taken to be the range in the parameter needed to make x 2 = Xmin+1> Xmin+4> 
and Xmi n + 9, respectively). Each of the parameters shown in the figure are well determined, since 
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the curves have smooth shapes with clear minima. The non-zero eccentricity is significant at about 
the 2cr level, as is the indication of non-synchronous rotation (e.g. the parameter = P ro t/-P rb is 
different than 1.0 by « 2a). On the other hand, models with no accretion disk are ruled out at the 
level of several a since the best-fitting model with no accretion disk has a x 2 of more than 210. 

As a check, extensive fits were made assuming different values of the X-ray eclipse width 0. 
Supplementary Table 1 shows the best-fitting values of the inclination, secondary star mass, black 
hole mass, and x 2 f° r various assumed values of 6. Over the whole range of 6 (40° — 50°), the 
inclination changes from 71.5° to 77.3°, and the black hole mass only changes from 15.97 M Q to 
15.35 M . 

Comparison of ellipsoidal model with additional light curves: Supplementary Figure [5] shows 
all of the available optical light curves of M33 X-7 phased on the X-ray ephemeris 5 : B and V 
light curves from the DIRECT survey (see ref. 10), g' and r' from the Gemini North Telescope, 
r' and %' from the Canada-France-Hawaii Telescope (CFHT) M33 Variability Survey (see ref. 11) 
and B, V, and I from WIYN. The CFHT and WIYN light curves have anomolously faint points 
near phase (I.e., near X-ray eclipse) that we cannot explain (there are no obvious defects in the 
images in question). These same light curves also have deeper minima near phase 0.5 (when the 
X-ray source is in front of the O-star companion), which, unlike the low points near phase 0, can 
easily be explained qualitatively by increasing the radius of the accretion disk from 45% of the 
black hole's Roche lobe radius to ~ 80%. The out-of-eclipse X-ray flux of M33 X-7 is variable, 
so an accretion disk with a changing radius is perhaps not unexpected. Alternatively, if the plane 
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of the accretion disk is not exactly parallel to the plane of the binary orbit, and the disk precesses, 
then its cross sectional area on the plane of the sky would change. This would result in changes in 
the depth of the light curve minimum near phase 0.5. 

We note that in the DIRECT, Gemini, and CFHT light curves the maximum near phase 0.75 
is consistently higher than the maximum near phase 0.25. The small orbital eccentricity in our 
model explains the slight difference between the heights of the maxima. 
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3 Supplementary figures, legends, and table 




Spatial position (arcsec) 



Supplementary Figure 1 | Spatial profiles of M33 X-7 and a close pair of neighbour 
stars. Shown are the profiles (filled circles) of the spectrum in the spatial direction near 
Hf3 of M33 X-7 (taller profile) and the neighbour stars for three observations with different 
seeing conditions: m 0.8 arcseconds full width half maximum (FWHM), left panel, « 0.6 
arcseconds FWHM, middle panel, and ss 0.4 arcseconds FWHM, right panel. Most of the 
observations were similar to the one shown in the middle panel. The solid line is a double 
Gaussian fit to the two profiles. The dashed lines show the individual Gaussians. In all 
cases, the overlap region between the two Gaussians is relatively small. The area in the 
overlap region is, from left to right, 1 8.6%, 6.3%, and 2.4% of the area under the profile of 
M33 X-7. The open circles show the profile of the residual image computed by SPECRES 
(e.g. the reconstructed spectrum minus the observed spectrum). 
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Supplementary Figure 2 | Comparison of GMOS pipeline and SPECRES extractions. 
{top) Typical one-dimensional spectrum as extracted using the GMOS pipeline (offset up- 
wards by 200 units for clarity) with nebular lines (e.g., H I, He I, and [O III]) present. 
{bottom) Spectrum produced using SPECRES. With the exception of H/3 and the stronger 
[O III] line, the nebular lines have been cleanly removed; note, however, the increased 
noise level. The regions of the spectrum falling in the detector gaps are indicated by the 
vertical lines. 
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Supplementary Figure 3 | Determination of the X-ray eclipse angle ©. The portion of the 
ACIS light curve from Figure 2 near the egress phase is shown with phase expressed in 
units of degrees. The horizontal bars represent the phase interval of the integration. The 
vertical error bars are la (s.d.) statistical. Different plotting symbols represent different 
ACIS Observation IDs: #6378 filled circles; #6382 filled squares; #7171 filled triangles; 
#7196 filled stars; and #7199 open circles. The thick solid line is the model used to 
determine 6 and is given by y(x) = c x for x < x 1 and y(x) = (c 2 - c l )/(x 2 - x x )x + c x for 
x > xi. The vertical lines denote, from left to right, eclipse widths of 40°, 46°, and 53°, 
respectively. For the full duration of the eclipse by the photosphere of the O-star, we adopt 
= 46°. 
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Supplementary Figure 4 | x 2 vs. parameter of interest. Projections of the \ 2 hypersur- 
face projected along various axes are shown (filled circles). The rr-axes give the value of 
either a fitted parameter (e.g. the inclination) or a derived parameter (e.g. the black hole 
mass), and the y-axes gives the optimal x 2 value when the corresponding rr-axis param- 
eter is fixed at its particular value. The dashed lines denote the formal la, la, and 3a 
confidence intervals. 
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Supplementary Figure 5 | Phased light curves and ellipsoidal models. All of the available 
optical light curves phased on the orbital period are shown with the ellipsoidal model fitted 
to the DIRECT B and V and the Gemini g' and r' light curves (solid lines). The dashed 
line shows an ellipsoidal model with a larger accretion disk that qualitatively explains the 
deeper minima near phase 0.5 in the CFHT r' and i' and the WIYN B, V, and / light 
curves compared to the DIRECT and Gemini counterparts. The error bars are la (s.d.) 
statistical. 
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(deg) 


% (deg) 


M 2 (M.) 


M(M ) 


x 2 


40.0 


71.50 


69.72 


15.97 


181.850 


41.0 


71.85 


69.70 


15.90 


181.132 


42.0 


72.36 


69.79 


15.79 


181.018 


43.0 


73.05 


69.89 


15.82 


180.460 


44.0 


73.51 


70.15 


15.77 


180.365 


45.0 


74.05 


70.08 


15.73 


179.914 


46.0 


74.62 


70.00 


15.65 


179.675 


47.0 


75.51 


70.24 


15.59 


179.637 


48.0 


75.84 


70.28 


15.49 


179.542 


49.0 


77.05 


70.13 


15.39 


179.546 


50.0 


77.27 


70.08 


15.35 


179.938 



Supplementary Table 1 Results for different values of X-ray eclipse width 6. The best- 
fitting values of the inclination i, mass of the secondary star M 2 , the black hole mass M, 
and the x 2 of the fit for various assumed values of the X-ray eclipse width are shown. 
The row in boldface is the one adopted in the main text. 
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